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ABSTRACT 

This paper discusses the efficient numerical treatment of optimal control problems gov- 
erned by elliptic PDE’s and systems of elliptic PDE’s, where the control is finite dimensional. 
Distributed control as well as boundary control cases are discussed. The main characteristic 
of the new methods is that they are designed to solve the full optimization problem directly, 
rather than accelerating a descent method by an efficient multigrid solver for the equations 
involved. The methods use the adjoint state in order to achieve efficient smoother and a 
robust coarsening strategy. The main idea is the treatment of the control variables on ap- 
propriate scales, i.e., control variables that correspond to smooth functions are solved for on 
coarse grids depending on the smoothness of these functions. Solution of the control prob- 
lems is achieved with the cost of solving the constraint equations about two to three times 
(by a multigrid solver). Numerical examples demonstrate the effectiveness of the method 
proposed in distributed control case, pointwise control and boundary control problems. 
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1 Introduction 


Computations in optimal control for distributed parameter systems demand huge com- 
putational resources. These are optimization problems with constraints which are par- 
tial differential equations. This paper focuses on the computational aspects of optimal 
control problems governed by elliptic partial differential equations or elliptic systems, 
where the control variables are of finite dimension. 

The classical optimal control algorithms use the costate ( a Lagrange multiplier) in 
order to define an iteration that converges to the minimum. These are descent type 
algorithms which involve the solution of the state and costate equations per iteration of 
updating the control variables. The state and costate equations, being discretizations 
of partial differential equations are quite expensive to solve. The over all solution of 
the optimization problem becomes therefore very expensive. 

The need for more efficient methods is obvious. One way of improving the situation 
is by using fast solvers , i.e., multigrid methods, for the state and costate equations at 
each step of a descent algorithm. This may result in a significant saving in computa- 
tional time, but will not affect the number of outer iterations for updating the control 
unknowns; it will just make the time required for each iteration smaller. 

A much more effective solution process can be obtained by having an algorithm 
that is aimed at the full optimization problem directly, rather than accelerating a 
linear equation solver. This paper is devoted to this approach. Multigrid algorithms 
that were designed to tackle the full problem rather than being a mere fast linear solver 
have been developed for a class of problems. The first one to be developed was the 
FAS algorithm for nonlinear equations [Bl]. Recently such an approach was successful 
in treating some bifurcation problems [T2] and stability calculation [T2], In all these 
cases maximal efficiency has been obtained as a result of treating the full problem 
directly. 

The high efficiency in the problems treated here is achieved by working on all 
unknowns of the problems (i.e., state, costate and control unknowns) simultaneously 
where scales of the different unknowns are taken into account. The main rule is to 
treat different control variables on different grids depending on their scale of influence. 
This means that a variable that has a non local effect on the solution should be treated 
mainly on coarse levels. On the other hand, a variable that has a non-smooth effect 
on the solution in some neighborhood, will be relaxed in that neighborhood on fine 
levels as well. More generally, local work on fine grids will be done to update a certain 
control variable only at the vicinity where it has a non smooth effect on the solution 
while work on coarser levels for that variable will be in larger and larger neighborhoods 
of that region. 

As is done in the classical optimal control algorithms a costate variable is introduced 
in our solution process. This allows the construction of good relaxation schemes in 
general, and is even more important when coarsening is considered (as explained in 
section 4). In the relaxation part of the algorithm we distinguish two main steps. One 
is designed to smooth the errors in the state and costate for a given choice of the control 
variables. This is done for the state variable using the Gauss-Seidel (GS) relaxations 
for scalar elliptic equations of state and Distributed-Gauss-Seidel (DGS) for systems of 
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elliptic equations [B2], The equations obtained for the costate are very similar to the 
state equations and are relaxed using essentially the same relaxation method as for the 
state equations. The other step of the relaxation is designed for updating the control 
variables. This step of the relaxation is done taking into account the smoothness of 
the change introduced in the state and costate as a result of a small change in one of 
the control variables. That is, some work is to be done on coarse grids only and some 
on fine levels in special subdomains. The update of the control variables is done in the 
following way. A small change is introduced in them, then its effect on the state and 
tie costate is computed approximately. The amplitude of that change in the control 
is then chosen so as to reduce the functional or to minimize the residuals of certain 
equations. 

The algorithms constructed that way are aimed at the full problem directly and 
need no iterations within iterations as the standard algorithms require. 

Numerical examples show that solutions to the levels of discretization errors are 
obtained in just a few work units, where a work unit is the work involved in relaxing 
the state equations. The complexity of solving a control problem is about twice that of 
solving the corresponding constraint equations, which is essentially optimal. Examples 
for distributed control, boundary control and pointwise control are given and all show 
the same efficiency. 

2 The Problem 

Let V, W,U and H be Hilbert spaces. Let A, B and C be operators as follows A : V — ♦ 

W, B : U — > VV, C : V — ► H and let U a d C U. Consider the problem 

nrin ue « td \\Cx — <f||^ . . 

Ax = Bu + f K } 

Here A is assumed to be an elliptic partial differential equation (PDE) or an elliptic 
system of PDE’s with smooth coefficients (see [ADN] ). B is assumed to be a finite 
dimensional operator, that is, there exist functions Tpj E W j = 1, , . . , q such that 

9 

B u — ^ min ( 2 . 2 ) 

i=i 

We assume that the functions ifij are smooth except for a finite number of singularities. 
U is then identified with lR q . The adjoint of f?, B* : W — ► IR q is given by 

B (f> — (<C Tpl j(f> >vv J • ' * > ^ ^ W ) (2*3) 

The following equations hold at the minimum, 

Ax = Bu + / 

A*p + C*Cx = C'd (2.4) 

B*p = 0 
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where B and B* are given by (2.2) and (2.3) respectively. Let xo be a solution of 
Ax o = / and d = d - Cx 0 . Problem (2.1 ) can be written in terms of x = x - x 0 and 
u as 


min ueWld ||C'i-d||J < (2.5) 

Ax = Bu 

So without loss of generality we can assume that / = 0. 

2.1 Descent Methods 

Iterative methods for solving the control problem are described. A triple (x,u,p) will 
be called compatible if it satisfies the equations 


Ax — Bu = 0 (2 6") 

A*p + C*Cx = C*d K ' 

Given an approximate solution (x,u) of (2.5), an improved one is obtained as follows. 
Denote by E(x,u ) the functional \\Cx - d|| 2 . For a compatible triple (x,u,p) 

E{x, u) = - < Cx, Cx > n -2 < u, B*p > u + <d,d> n (2.7) 

We look for a change 71Z in u and corresponding changes in x,p that result in a new 
compatible triple, and for which the value of E is smaller. Let (2 -f* 7®> ® d - P d" TP) 
be a new compatible triple. Then 


E[x + 72,11 + 7 u) = E(x , it)+ 

-27 < Cx,Cx > H -7 2 < Cx,Cx > n (2.8) 

-27 < B*p,u > u -27 < B*p,u > u -2-y 2 < B*p, u > n 


Choosing u = B*p , and x,p satisfying Ax — Bu — 0 and A*pA C*Cx — 0 we get 


E(x + 71,11 + yu) = E(x, u) 

— 2 7 (< Cx,Cx >„ + < B*p,B*p > u + < B*p,u > u ) (2.9) 

-7 2 (2 < B*p,B*p > u + < Cx,Cx > H ) 

Among all possible 7, the choice 


1 < Cx,Cx > n + < B*p,B*p > u + < B*p,u > u , 0 1fV) 

7 _ "2 \\Cx\\l+2<B*p,B*p> u 

minimizes the functional in the direction chosen. This can be a basis for an algorithm 
for solving the optimal control problem. 


A Descent Algorithm 
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Set uo = 0 ,po = 0,7 = 1 , n = l, 

While 7 > e Do 
Begin 

1 . Solve for x n the equation Ax n — Bu n = 0 

2. Solve for p n the equation A*p n + C*(Cx n — d) = 0 

3 . Set u = B*p n . Solve Ax — 55 = 0 , yl*p + C*Cx = 0 . 

(approximately) for 4 ,p respectively. 

4 . Set u n + 1 = u n +7u, x n+i = x n + 72, p n+1 = p n + 7 P» n = n + 1 
with 7 given by (2.10), 

End 

Approximate Descent Algorithm 

Set uo = 0 ,po = 0,7 = l,7i=l. 

While 7 > e Do 
Begin 

1 . Starting with x n = x n _i relax Ax n = yielding 

a solution satisfying ylx n — Bv^ = 

2. Starting with p n = p n -i relax A*p n + C*Cx n = C*d, 
yielding a solution satisfying A*p n + C*(Cx n — d) = r” 

3 . Set u = B*p n . Solve Ax — Bu — 0 , A*p + C*Cx = 0 
(approximately) for x,p respectively. 

4 . Set u n _|_i = u n 7U, Xn+i = x n -f 7^ » Pn- {-1 = Pn 4 * 7 P* u = 71 -(- 1 
with 7 given by (2.10). 

End 

e is a small enough number, depending on the accuracy required for the solution. 
Note that if ||r£|| — > 0 and ||r-”|j -* 0 then the above algorithm converges to the 
minimum of ( 2 , 5 ). 

The expensive steps in this algorithm are the ones for obtaining x and p by ap- 
proximately inverting A and A* } which are elliptic operators (or systems). Moreover, 
the solution of the elliptic equations governing the state and costate needs to be done 
many times. 

Remark Considering the form of the necessary conditions one may be tempted 
to choose 7 based on minimizing the L 2 residuals of the third equations of ( 2 . 4 ). This 
leads to 


< B*p,B*p > 

< B*p,B*p> 


( 2 . 11 ) 
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which may not be a good choice as < B*p,B*p > may vanish, causing the iteration to 
get stuck away from the minimum. 


3 Naive Multigrid Approaches 

The approximate-descent algorithm presented in the previous section suggests an im- 
mediate acceleration procedure, namely, using some fast solver in steps 1,2, and 3. 
Since the form of the operator A* is very similar to that of A the same fast solver can 
be used in all the steps. The operator A is assumed to be elliptic and therefore calls 
for multigrid methods. 

Assume that we are given V k , Wfc , U k and 7i k Hilbert spaces approximation the 
original V,W,U and H. Let A k ,B k and C k be operators as follows A k : V k -+ W k , 
B k -Uk -* W*, C k : V k -> H k . Assume also that interpolation and restriction operators 
P%,P£,Pl;,Rx>Rp> R t are given; R £ : V k +i -» V k , Rp ■ VVfc+i -» W fc , R * : Z4+i -*• Uk, 
P% : V fc -» Vjb +l , Pp : Wfc -» Wfc+i, P$ : U k -* U k+1 . The superscript k in these 
expressions represent a level of discretization where k — 1 correspond to the coarsest 
level. 

Denote by MG(A k , b, x, y, i/ a , i/ 2 ) a multigrid V(vi, v 2 )-cyc\e for solving A k x = b on 
level k, starting with initial approximation x ending with y. Assume that a solution of 
the minimization problem is needed for level k = m. A possible fast algorithm for it, 
starting with initial approximation (x m ,u m ,p m ) is denoted by 

(x m ,u m ,p m ) <- MinMG(A m ,B m ,C m ,x m ,u m ,p m ,cr ) (3.12) 

and given by 

Set = 0 ,P o m = 0 , 7 =l,n = l. 

While 7 > e Do 
Begin 

1 . Perform a cycle MG(A m ,B m u™,x™,x™ + 1 ,vi,V2) 

2 . Perform a cycle MG(A^ nt C^ n d m - C^Cx™ ,p™ ,p™ +i , i'i, ^2) 

3 . Perform a cycle 5 m 5 ^y m , i™ , 1™ ±\,v\,v 2 ) 

4 . Perform a cycle MG{An,-CmCm*™+\,i% 1,^2) 

5 . Set < +1 = < + 7^+1 . C+i = ®n + TC+1 . Pn+1 = Pn + 7 #h-i • 
n = 7 i + 1 and 7 given by (2 . 10) . 

End 

An obvious improvement of this algorithm is to start with a good initial approxi- 
mation for the optimization problem on level m. This can be achieved easily by solving 
the problem first on level 771 — 1. Applying this idea recursively one arrives at the 
following algorithm. 
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Set Uq — 0,Xq = 0,Pq = 0 ,k ~ 1 . 

While k < m Do 
Begin 

1. Perform a cycle 

(: x k ,u k ,p k ) ♦- MinMG(A k ,B k) C k ,x k ,u k ,p k ,d k ) 

2. If k = m stop 

3. x k+1 = E k x k ,p k+1 = n k p k ,u k+1 =Piu k ,k = k + 1. 

End 

where II*, lip are operators analogous to P* and P* , respectively, but of possibly a 
higher order (as is usual in multigrid algorithms for elliptic problems)* 

Our aim is to develop an algorithm in which the total computational cost of solving 
the control problem will be only few times (2-3) more than that of solving for the state 
variable x given that u is known. This is done in the next sections where a multigrid 
approach which is aimed at the full problem directly is constructed. 

4 “One Shot” Multigrid M eth ods 

In this section we describe a multigrid method that is aimed at the full optimization 
problem, rather than serving as a fast solver in a step of a basic optimization algorithm, 
as was done in the previous section. 

The construction of this method has been achieved basically by following the next 
five steps. 

1. Distinction between smooth versus non-smooth components in all types of vari- 
ables involved 

2. Construction of a basic non-expensive relaxation for the full problem. 

3. Classification of rate of convergence of the basic relaxation for the different vari- 
ables in the problem, with regard to their smoothness. 

4. Based on 3, determination of the role of coarse grids in accelerating the fine grid 

convergence. : : ^ : 

5. Construction of a coarse grid problem that approximate the Full fine grid residual 
problem. 

Steps 1,2 and 3 are discussed in detail in section 4.2, the rest in section 4.1. The 
resulting algorithm involves a sequence of discrete optimization problems, starting 
from a very coarse discretization ending with a fine grid one on which a solution of the 
problem is required. The optimization problem on each level is served to accelerate the 
convergence of the next finer level problem. The iterative process for solving the finest 
grid problem involves a relaxation process on each level whose effect is to damp out the 
non-smooth part of the error, i.e., the part of the error which cannot be represented 
on the next coarse grid in the sequence. This together with a proper transferring of 



information (residuals and current fine grid solution) from fine to coarse levels and 
back to fine levels results in a very high efficiency. A full description of the method 
includes a relaxation scheme and a coarsening strategy. 


4.1 Coarsening 

A coarsening scheme requires the definition of coarse grid spaces analogous to V, W, U, U 
which we denote by the same letters with a superscript c. Also inter grid transfer 
p x> p p> p u and R x , Rp, Ru and coarse grid operators A c , B c , C c are required. Having all 
this we can define the coarse grid minimization problem. An attempt to use 


minuc^ \\C f {x f + P x x c ) - d f \\ 1 2 ( 4 . 1 ) 

A c x c - B c u c = R x (d ~ A f x * + B f uf ) 

as an approximation for the error on the fine grid fails. The reason for this is that 
this coarse grid minimization problem does not preserve the fine grid minimum once 
it was reached. This follows immediately from the necessary conditions of the coarse 
grid problem (4.1), namely, 

A c x c - B c u c = R x {g f x - A fX f + B f uf) 

A* c p c + C* c C c x c = R p C* f {df - C f xf) (4.2) 

B* c tf = 0. 

The coarse grid minimization problem formulated above is a problem for the correc- 
tion in the fine grid approximation. Therefore, upon reaching the fine grid minimum 
the coarse grid correction x c , u c should be zero. This however is not the case since the 
right hand side of the second equation is non-zero in general at the minimum of the 

fine grid problem. ....... 

In order to define a correct coarse grid problem we rewrite the fine grid minimization 

problem as 

min u/ u s < CfX f } x f > -2 < x*,gf, > -2 < u f >g£ > (4 3 ) 

Afxf - BfU f = gl 

where g{ = <7, 5 / = C*d and gl = 0 on the fine grid. Note that the problem is defined 
once g£,gl, g f p are given . 

1. Correction Scheme. The coarse grid problem for the correction is of the same 

form, that is, 


nun 


:eM . £ d 


< C c x c ,x c > -2 < x c ,Pp > -2 < it c ,<7u > 


A c x c 


B c u c = g x 


(4.4) 
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where 


9 c x = R x (gl-A f x' + B f uf) 

9% = »p(sl ~ AjP f ~ C* c C c xf) (4.5) 

9 C u = RM ~ B}pf). 

The necessary conditions for this coarse grid problem are 


A c x c - B c u c = g% 

A* C P C + C*C c x c = g‘ (4. 6 ) 

KP C = 9 C u 

which are the Correction Scheme (CS) (see [B2]) equations for the fine grid necessary 
conditions. In particular they have the property of introducing no change to the fine 
grid problem once the minimum has been reached. 

Once an approximation to the coarse grid problem is found it is used to correct the 
fine grid approximation by 


<— x* + P x x c 

u f <- u f + P u u c ( 4 . 7 ) 

P* *- P f + PpP c 

Note that this coarse grid problem depends on the costate p obtained after the fine 
grid relaxation process. Without using it the coarsening becomes very difficult as was 
explained before. 


2. Full Approximation Scheme. The coarse grid problem for quantities that 
approximate the full fine grid solution rather than the correction for it have similar 
form. The minimization problem is 


< C c x c ,x c > -2<x c , g ;> -2< u c ,g c u > 

A c x c - B c u c = g% 

where 

9% = Rx{g f x - Afx' + Bfuf) + A c R x xf - BcRvuf 
9p = Rp{g f p - A)pl - C*C c xf) + A* c Rppf + C* c C c R x xf (4.9) 

9t = Ru{gL ~ B)p f ) + 5 c %p/. 

The operators R x , R p , are fine to coarse grid transfers possibly different from R x , Rp, R u 
(see [B2]). The necessary conditions for this coarse grid problem are 


A c x c - B c u e = g% 

Atf + CtC^^gi 

B' c P c = gl 


(4.10) 
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which are the FAS equations for the fine grid necessary conditions. In particular they 
have the property of introducing no change to the fine grid problem once the minimum 
has been reached. 

Once an approximation to the coarse grid problem is found it is used to correct the 
fine grid approximation by 


x* *- x i + P x {x c - R x x f ) 

uf <- uf + P u (u c - Rvul) (4.11) 

pf *-p f + P p {p c - RpP f ) 

4.2 Relaxation 

The general equations for state and costate variables on any of the grids involved in 
the process have the form. 


Ax — Bu — f x + 

A*p + C*Cx = f p [ * ' 

These two equations are being relaxed on all levels using one of the standard relax- 
ation methods for elliptic problems. That is, Gauss-Seidel or damped Jacobi in case A 
is a discretization of a scalar PDE, and Collective Gauss-Seidel (CGS) or Distributed 
Gauss-Seidel (DGS) when systems of elliptic PDE’s are involved. So the problem of 
relaxation for x and p is quite standard. An important property of such relaxation 
schemes is that for a wide class of discretization methods ( h-elliptic discretizations) 
they have fast convergence for h-non-smooth errors (defined later) (see [B2]). This 
implies that the errors in the x and p equations will have smooth residuals after just a 
few relaxation sweeps. This does not mean that the errors for the optimization prob- 
lem are smooth. The relaxation of equations (4.12) serves as a process that brings the 
solution closer to the constraint surface defined by the first equation there. 

As far as updating the control variables the situation is more complicated. Note 
that a small change in any of the control unknowns Uj implies a change in x and p 
everywhere in the domain. Good efficiency for the full algorithm is to be obtained if the 
changes in x and p that result from changes in u, can be calculated with a computational 
cost which is not much greater than that of relaxing the equations (4.12) once. 

The smoothness assumption on the coefficients of the operator A implies that non- 
smooth errors correspond to non-smooth residuals and vice versa. Non-smooth errors 
are also fast to converge, for h-elliptic discretizations, as mentioned before (this is not 
the case for quasi-elliptic discretization schemes). This results in the following basic 
rule for updating the control variables. 

• MaJfce changes in the control variables that result in non-smooth residuals 

Note that a change in any of the control variables may result in smooth changes 
on one grid but non-smooth changes on a coarser grid. This suggest that the different 
control variables will be relaxed on possibly different grids depending on their scale of 
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influence in the solution. To quantify more precisely these ideas we need the following 
definitions. 

Definition: A function V'/iCO twM be called h-non-smooth at £ £ iff 

\(phi>h)(0 - > a#h(£)l o < a 0 < a < l,a 0 = 0(1), (4.13) 


where ph is a local averaging operator. 

Definition: A function iph(£) wz7Z be called uniformly h-non-smooth in fl iff it is h- 
non-smooth in every point £ 6 Q h 

Thus if ipj is uniformly h-non-smooth its effect on the solution x and p can be 
calculated accurately with just a few relaxation sweeps. Typically 2-3 relaxation sweeps 
are enough to reduce the errors by an order of magnitude. This will be the key to our 
method of updating the control variables Uj. Let 

Hy = £ Q h : is h-non-smooth at (4-14) 

If h correspond to the coarsest level and Hy / p then we define Hy = Define 


D k (() = diagtif^f) *,(()) 


< 4.(0 


= {o < 


if Tp% is h-non-smooth at £ 
otherwise 


(4.15) 


The qxq matrix D h will be used in defining the relaxation for the control variables. 

In case the functions ^y are uniformly h-non-smooth the update of the control 
variables becomes much simpler. The perturbation in u that is done in minimizing 
the functional is by 7 D h (B h )*p h , followed by a few relaxation sweeps of the state and 
costate equation in a vicinity of IJEiy, yielding a good enough approximation to x h 
and The actual perturbation is then chosen so as to minimize the functional or the 
other choice as explained in section 2. Note that the multiplication by D ^ causes the 
perturbation in the right hand side to be h-non-smooth on the given level. That is, 
each control variable is being updated on the appropriate grid. 

For ipj with a singularity a slightly different approach is taken. Here the effect on the 
solution is smooth away from the vicinity of non-smoothness of ^y. This follows form 
basic theories of elliptic partial differential equations [ADN]. Thus, if the observation 
operator C is supported far enough from the singularities of $y, even coarse grids can 
compute accurately the effect on Cx, which is the important quantity in updating 
the Uj. However, when the singularities are close to the support of C this is not the 
case and some refinement has to be done locally in order to account for the local non- 
smoothness. In practice, a few points around the singularity need to be relaxed, further 
points are relaxed on coarser grids. The jfinest grid which needs to be involved in the 
process with a given singularity is determined from the accuracy achieved in Cx which 
can be estimated using the quantities 

|| C h x h - C h x 2h || (4.16) 


and the knowledge about the accuracy of the discretization involved. The later, if not 
known, can be computed from 


jjcW - C h x 2h || 
|| C h x 2h - C h i 4,l || 


(4.17) 
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In any case, the work performed in the calculation of should not involve more 

than one multigrid cycle, done as an FMG cycle used with local refinements. The exact 
formulation of the algorithm for updating the control variables in the most general case 
will be presented elsewhere. We denote the process of updating the control variables 
in the direction u by 

(x,p,u) <— MinRe^XjP^u^u) (4*18) 

Having defined the relaxation and the coarsening process the full algorithm can 
now be described. 

4.3 Unistep Algorithms 

Let Vk,WkMk and H k be Hilbert spaces approximating V, W,W and H , respectively. 
Let Ak y Bk and Ck be operators as follows A k : V* — ► W*, B k :U k -* W*, C k : V* — > 
Tiky approximating the operators A,B,C, respectively.^ Assume also that interpolation 
and restriction operators P k ,P k ,P k ,R k x , R k , R k , R k x ,R$, are given as R k : V k +x -*■ 
V k , R k : Wjfe+i -► W fc , K ■ *4+i -» *4, Rl ■ Vk+i -* V ft> R$ : Wa=+i -» W fc> 
St : £4+1 - Uk, P x fc ■ Vk - V fc + 1 , P k : W fe -► Wjt+i, P u fc : *4 - Wfc+i- Also assume the 
the qxq matrices D* (defined in section 4.2) are given. 

On all levels minimization problems as follows are given, 

min u fc eW i. d < Ck x k ,x k > -2 < x k ,g% > — 2 < u k ,gt > (4-19) 

A k x k = B k u k + g k 

Observe that this is equivalent to having the following systems of equations on all 
levels 


A k x k = B k u k + g k 

A* k p k + C* k C k x k = 5 p ( 4 - 2 °) 

B* kP k = gt ■ 

Let ( x k ,p k ,u k ) be an approximate solution of (4.19). We define next an algorithm 
for improving it, denoted by 


(x k ,p k , u k ) 4 - MG(x k ,p k ,u k ,g k x ,g$,gt) 


(4.21) 


If As = 1 Then 

1. relax the first two equation in (4.20) until convergence. 

2. iterate until convergence: 

2a. perform the cycle 

{x k ,p k ,u k ) 4 - MinRel(x k y,u k } D k (gt-BkP k )) 
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2 b. relax the first two equations in ( 4 . 20 ). 

Else 

1 . perform the following v\ times 

la. relax the first two equation in ( 4 . 20 ) 

lb. perform the cycle 

{x k ,p k ,u k ) «- MinRel{i k , 1 f i u k ,D k {g k -B' k p k )) 

2. Let k = k — 1 , and 

9x = R* (^ +1 - A k+1 x k+1 + B k+1 u k+1 ) +A k R k x k + 1 
9p = R k ($ +1 ~ At +1 p k +' + C* k+1 C k+1 u k +') + A k R$p k +' 
ffu = R k P (s* +1 - B£+iP k+1 ) + B k R$p k +i 

3. perform 7 times the cycle 

{x k ,p k ,u k ) *- MG{x k ,p k ,u k ,f k ) f k ,f k ) 

4. correct fine grid solutions 

z* +1 <_ x k+i + _ p* j.fc+1^ 

p fc+1 <- p fe+1 + P k (p k - R k p k+1 ) 

u k+1 u *+i + _ R k u k+ij 

5. perform the following v\ times 

5a. relax the first two equation in ( 4 . 20 ) 

5b. perform the cycle 

(x k ,p k ,u k ) <- MinRel(x k ,p/° ,u k , D k (g k — B k p k )) 

End 

In order to obtain full efficiency the algorithm starts at the coarsest level, where 
each time a refinement is done followed by a fixed number of MG cycles. This is the 
N-FMG algorithm which is defined next. 

1. Solve (4.19) for k = 1 

2. It = k + 1, x k = II*- 1 !*- 1 , p k = II k ~ l p k ~ x ,u k = P k - l u k ~' 1 . 

3. Define g k ,g k ,g k . 

4. Perform N time the cycle 
(x k y,u k ) «- MG{x k ,p k ,u k ,f k J$,f k ) 

5. If k = M stop, else goto 2 . 

At the end of this algorithm an approximate solution to the minimization problem 
on level k = M is given. 

5 Numerical Examples 

Numerical experiments were conducted with scalar elliptic problems governed by the 
Laplacian in the unit square with Dirichlet boundary conditions. That is, using the 
notation of section 2, 
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(5.1) 


ft = {(i, y) : 0 < x,y < 1} 

V = W = {<£ e £ a (ft) : <t>\an = 0} 

A = A 

h = l 2 ( r) r = {(*,o) :0<x< 1 } 

C = d/dn\r 

The control space was M q . Different cases were considered for (ipj j — 1 , . . . , g) , 
in the definition of B. 

Uniform grids were used in the discretization with the standard 5-point formula for 
the Laplacian, namely, 


ft* = {(i/N,j/N) : 0 < i,j < N} 

T h = {(i/N,0):0<i<N} (5.2) 

(aW)u = M x Ui + x M.i + x t,-i + *&+i - «&) 

where Nh = 1, h being the discretization parameter. The observation operator was 
discretized using a 3-point formula, 

(C"**' 1 )^ = i(1.5JT{b - 2X>i + .5X& (5.3) 


and the functions ^ were discretized by simple injection. 

With this discretization the discrete solution is expected to have second order ac- 
curacy. That is, 

h -X\\ <(3 x h 2 
h -u\\< /3 u h 2 

where /3 X ,/3 U are constant independent of h (they depend on high order derivatives of 
the state variable). It is enough to solve the discrete problem to an accuracy satisfying 




||* fc -x fc ||<pr fc -x|| (55) 

||v* _ u *|j < || u * _ tt || k ; 

Here the tilde quantities represent the current numerical approximations to the solution 
of the discretized problem. The quantities ||X* — X||, \\u h — ii|| are called discretization 
errors. 

In all the examples reported below the algorithm described in section 5 was used 
with the following parameters: v\ = 2,u 2 = 1. The restriction operators for both the 
residuals and the full solutions R^R*, Rx,Rp were the 9-point full weighting, that is 
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1 2 1 

2 4 2 
1 2 1 
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(5.6) 



The interpolation operators P* +1 , P* +1 were the standard bi-linear interpolation, II*, Il£ 
the bi-cubic interpolation and R* = R u = I , the identity operator. The mesh size on 
level k was 2~(* +1 ), k=l being the coarsest level. Results for 3-FMG-V(2,l) cycle are 
given. In each table results for the L 2 residuals r x ,r p ,T u of the state, costate and the 
control equations are given. Their exact definition is given by 

Nk 

M2 = £ ((ai)ij + (B t u% - (A t x%fhl (5.7) 

i,j=l 

IMI2 = E ((«*)« - (5.8) 

»,i=i 

Ml* = E((sj)l - (Biv k )>)* (5-9) 

i=l 

where the scaling of h ^ in these definitions is used so that these norms approximate 
the continuous norms and so residuals on different levels can be compared. 

1, Distributed Control 

In the distributed control case we have worked with the following problem 


min f (tt- — d) 2 do 
« Jr on 

(5.10) 

Ax = £i=i Uii>i 


x|an = 0. 

(5.11) 


ipj - sin(nj'K^i)sin(mjT^ 2 ) n j, m j integers, 

Table 1 shows the results for a on e- dimensional control case, that is, q = 1 where 
m\ — nj = 1, i.e. a smooth case. This is the simplest of all cases and is given here 
mainly for reference. Residuals of the state, costate and control equations are given 
separately. The error in the control (relative to the true differential error) is given as 
well. Observe that a 1-FMG-V cycle gives solutions to the levels of discretization errors 
on fine levels, as ||u - u eiact || reaches its minimum in essentially one cycle. The 0(h 2 ) 
behavior of the error is clear from the results, reflecting the order of the scheme used. 

Table 2 shows the results of a similar experiment in which the control is three 
dimensional (q = 3) with n\ = n 2 = 1, w,\ = 7712 = l,n 2 = m 3 = 2. The behavior here 
is similar to the previous case. The 0(h 2 ) is evident also here. 

2. Boundary Control 

The boundary control case problem was 

min f - d) 2 do (5.12) 

« Jr on 
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(5.13) 


Ax — 0 

^len = E?= i 

ipj = sin(nj7r£i)|an nj integer, 

Tables 3 and 4 shows results for this problem. The first with q = l,ni = 1, the 
other with q = 2,ni = l,n 2 = 2. Also in this case the results are basically optimal, a 
solution to the levels of discretization errors are obtained in 1-FMG-V(2,1). The 0(h 2 ) 
convergence toward the differential solution is evident in both cases. Observe that the 
initial residuals on finer levels start are smaller by a factor which is close to four. This 
is typical to examples with smooth solution and proper FMG interpolation operators 

n* , iip . 

3. Pointwise Control 

The pointwise control case problem was 


min f (— d) 2 da 

« Jr on 

(5.14) 

As = E’= 1 

(5.15) 

&|en = 0- 


tpj — ^ 

Tables 5,6 and 7 show results for non-smooth control. Results are given for one 
two and three delta functions, i.e., q = 1,2,3. The location of the delta function is 
given in each table. Here the solution is less smooth than before. Still the results 
show essentially the same behavior as in the smooth case. Observe that the initial 
residuals on the different levels are of the same order, reflecting the non- smoothness of 
the solution. This behavior cannot be improved by using a higher interpolation in the 
refinement stage of the FMG algorithm, although local relaxation in the vicinity of the 
singularities can improve the results we have not experimented with such ideas here. 

In some of the experiments it can be observed that one of the three residuals shown 
is increased in some of the cycles. This is reasonable since the three residuals are 
coupled and it is only their sum which goes down in each cycle. The results presented 
by all tables clearly demonstrate the effectiveness of the method developed here which 
aimed at the full optimization problem, therefore leading to one-shot solution of these 
problems. 1-FMG-V cycle is basically enough to reach below the levels of discretization 
errors. 
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level 

cycle no. 

■ 

rxh 

■■ 

IWU 

||u — U eX act | ( 2 

1 

5 

.337e— 06 

.121e— 05 

.397e— 09 

.8166 

2 

1 

.591e— 02 

,121e+00 

.415e— 03 

.589e— 01 


2 

.104e— 03 

.609e— 02 

.975e— 04 

.612e-01 


3 

.216e— 03 

.373e— 02 

.485e— 04 

.568e-01 


4 

.592e— 04 

.996e— 03 

.183e— 04 

.581e-01 


5 

.250e— 04 

.446e— 03 

.737e— 05 

.576e-01 

3 

1 

.860e-03 

,522e— 01 

.319e— 04 

.174e— 01 


2 

.648e— 03 

.292e— 02 

.209e— 06 

.149e— 01 


3 

.121e-04 

.136e— 03 

.755e— 07 

.147e-01 

4 

1 

.125e— 03 

.357e— 01 

.434e— 05 

.329e-02 


2 

.324e— 04 

.132e— 02 

.371e— 06 

.305e— 02 


3 

.407e— 04 

.683e— 04 

.294e— 07 

.301e— 02 


Table 1: q = 1 ,n x = l,m 1 — 1 


level 

cycle no. 

■HIM 

■ 

S19H 

BB 

^^2 

91 

1 

5 

.246e— 06 

.110e-05 

.374e— 08 

.253e+00 

2 

1 

.887e— 02 

.207e+00 

.654e— 03 

,174e+00 


2 

.378e-02 

.592e— 01 

.200e— 03 

.117e+00 


3 

.169e— 02 

.258e— 01 

t- 4 

O 

00 

C* 

1 

O 

CO 

,994e— 01 


4 

.237e— 03 

.359e— 02 

.271e— 04 

.lOOe+OO 


5 

.905e— 04 

.130e— 02 

.908e— 05 

.102e+00 

3 

1 

.246e— 02 

.llOe+OO 

.772e— 04 



2 

.418e— 03 

114e— 01 

.222e— 04 

.317e— 01 


3 

.113e— 03 

.294e— 02 

.755e-05 

■mi 

4 

1 

.380e-03 

.618e-01 

.835e-05 

.858e-02 


2 

.595e— 04 

.316e— 02 

.110e-05 

.630e— 02 


3 

.383e-04 

.167e— 03 

.166e— 06 

.618e-02 


Table 2: g = 3, 7 ii = n 3 = 1 , mi = m 2 = 1, n 2 = m 3 = 2 
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level 

1 

2 


3 


4 


cycle no. 
5“ 
1 
2 

3 

4 

5 
1 
2 
3 
1 
2 
3 


IMI2 
.631e— 06 
.466e— 01 
.159e-01 
.521e-02 
.177e— 02 
.597e-03 
,826e-02 
.184e-03 
.668e— 04 
.158e— 02 
.338e— 03 
.273e— 03 



.566e— 07 


.122e+01 
.379e+00 
.126e+00 
.426e— 01 
.144e— 01 
•480e+00 
.146e— 01 
.641e— 03 
.264e+00 
.199e— 01 
.253e— 02 


.000e+00 
,126e— 02 
.123e— 03 
,142e— 04 
,161e— 05 
.183e— 06 
.431e— 07 
.332e— 10 
.200e— 10 
.689e-06 
.285e— 08 
,895e-ll 


^ Inexact 1 2 

.179e+00 
.212e— 01 
.617e-01 
.477e— 01 
.524e— 01 
.509e— 01 
.139e— 01 
.133e— 01 
.133e— 01 
.225e— 02 
.279e— 02 
.275e-02 


Table 3: g = 1,71* = 1 


level 

cycle no. 

mmm 

mmm 

■HIM 

||w - u exact \\ 2 

1 

5 

.631e-06 

.800e— 06 

.296e— 13 

.334e+00 

2 

1 

.354e+00 

.786e+01 

.288e— 01 

,149e— 01 


2 

.206e+00 

.388e+01 

.896e— 02 

.150e+00 


3 

.120e+00 

.229e+01 

.307e— 02 

.754e— 01 


4 

0 

01 
ct> 

1 

0 

.134e+01 

.106e— 02 

.117e+00 


5 

.415e— 01 

.792e+00 

.367e— 03 

.930e— 01 

3 

1 

,523e— 01 

.250e+01 

.924e— 03 

.392e— 01 


2 

.106e-01 

.526e+00 

.446e— 04 

.298e-01 


3 

.212e-01 

.104e+00 

.168e— 05 

.281e-01 

4 

1 

.120e-01 

.126e+01 

.394e— 05 

.528e— 02 


2 

.628e-03 

.668e— 01 

.958e— 08 

.600e— 02 


3 

.274e-03 

.466e— 02 

.313e— 10 

.594e-02 


Table 4: q = 2,ni = 1 ,n 2 = 2. 
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level 

1 

2 


3 


4 


cycle no. 

5 

1 

2 

3 

4 

5 
1 
2 
3 
1 
2 
3 


IMa 

.477e-06 
.485e+00 
.395e— 01 
.320e-01 
.931e-02 
.364e— 02 
.341e+00 
.861e— 02 
.248e— 02 
.574e+00 
.132e-01 
.587e-03 


IIMa 
.853e— 07 
.172e+01 
• 167e+00 
.310c-01 
.131c-01 
.445e— 02 
.791e-01 
.277e— 01 
•510e-02 
.332e— 01 
,112e— 02 
.200e-03 


lku||a 

.279e— 08 
.363e— 02 
.332e-03 
.132e-03 
.447e— 04 
.166e— 04 
.115e-03 
.118e— 04 
.258e-05 
.364e— 05 
.770e— 07 
.883e— 07 


l|tt ~ Uexact 

.370e+00 
.536e-01 
.863e-01 
.572e— 01 
.656e— 01 
.624e— 01 
.114e-01 
.128e-01 
.138e-01 
.208e— 02 
.272e— 02 
.275e-02 



Table 5: q = 1,^ = (.5, .5) 


level 

~T~ 

2 


3 


4 


cycle no. 
_ 

1 

2 

3 

4 

5 
1 
2 
3 
1 
2 
3 


IMa 

.211e-05 
.619e+00 
.236e+00 
• 147e+00 
.927e— 01 
.598e— 01 
.391e+00 
.389e-01 
,138e-01 
.660e+00 
.168e-01 
.708e— 03 


\Kh 

.921e-05 
•437e+01 
.223e+01 
.144e+01 
.921e+00 
Jttle+OO 
.389e+01 
.200e+00 
.336e— 01 
361e+00 
.187e-01 
.120e-02 



.347e— 06 


.363e— 01 
.256e-01 
.166e— 01 
.105e— 01 
.677e— 02 
.148e— 01 
.160e— 02 
.523e— 03 
.587e— 03 
.405e— 04 
.116e— 06 


lh ~ u exact \\ 
.864e+00 
.390e+00 
.181e+00 
.536e— 01 
.425e— 01 
.909e— 01 
.537e-01 
.373e— 01 
.312e-01 
.673e— 02 
.511e-02 
.509e-02 


2 


Table 6: q = 2,6 = (.75, .5), 6 = (.5, .5) 
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BSSI 

cycle no. 

■IBM 

— 

— 

BUS HI 

i 

5 

.131e— 05 

.953e— 05 

.334e— 05 

.115c+01 

2 

1 

2 

3 

4 

5 

.816e+00 

.293e+00 

.215e+00 

•131e+00 

.110e+00 

.447e+01 

.245e+01 

.177e+01 

.119e+01 

.883e+00 

.477e— 01 
.319e— 01 
.243e— 02 
.162e— 01 
,131e— 01 

.531e+00 
.273e+00 
.848e— 01 
.345e-01 
.127e+00 


1 

2 

3 

.396e+00 

.431e-01 

.244e-01 

.412e+01 
.214e+00 
.476e— 01 

.155e— 01 
.236e— 02 
.103e— 02 

.793e-01 
.640e-01 
.531e— 01 

4 

1 

2 

3 

.683e+00 
,402e— 01 
,452e— 02 

.373e+00 
.341e— 01 
.533e— 02 

.132e— 02 
•258e-03 
.3096-04 

.168e-01 
.880e— 02 
.788e— 02 


Table 7: q = 3,6 = (-75, .5), 6 = (.5, .5)6 = (-25, .5) 


20 






















NASA 

Natonai Aewautcs anc 
Space Aorn*msif3lon 


Report Documentation Page 


1. Report No. 

NASA CR- 187497 
ICASE Report No. 


91-2 


2. Government Accession No. 


4. Title and Subtitle 

"ONE SHOT" METHODS FOR OPTIMAL CONTROL OF DISTRIBUTED 
PARAMETER SYSTEMS I: FINITE DIMENSIONAL CONTROL 


7. Authorls) 

Shlomo Ta T asan 


9. Performing Organization Name and Address 

Institute for Computer Applications in Science 
and Engineering 

Mail Stop 132C, NASA Langley Research Center 
Hampton, VA 23665-5225 


12. Sponsoring Agency Name and Address 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23665-5225 


15. Supplementary Notes 

Langley Technical Monitor: 
Richard W. Barnwell 

Final Report 


3. Recipient's Catalog No. 


5. Report Date 
January 1991 


6. Performing Organization Code 


8. Performing Organization Report No. 

91-2 


10. Work Unit No. 


505-90-21-01 


11. Contract or Grant No. 

NAS1-18605 


13. Type of Report and Period Covered 
Contractor Report 


14. Sponsoring Agency Code 


Submitted to SIAM Journal 
on Optimization 


16. Abstract 

This paper discusses the efficient numerical treatment of optimal control prob- 
lems governed by elliptic PDE's and systems of elliptic PDE T s, where the control is 
finite dimensional. Distributed control as well as boundary control cases are dis- 
cussed. The main characteristic of the new methods is that they are designed to 
solve the full optimization problem directly, rather than accelerating a descent 
method by an efficient multigrid solver for the equations involved. The methods use 
the adjoint state In order to achieve efficient smoother and a robust coarsening 
strategy. The main idea is the treatment of the control variables on appropriate 
scales, i.e., control variables that correspond to smooth functions are solved for 
on coarse grids depending on the smoothness of these functions. Solution of the 
control problems Is achieved with the cost of solving the constraint equations about 
two to three times (by a multigrid solver) . Numerical examples demonstrate the ef- 
fectiveness of the method proposed in distributed control case, pointwise control 
and boundary control problems. 


17. Key Words (Suggested by Authorls)! 

unistep methods, distributed parameter 
systems, multigrid, control 


18. Distribution Statement 
64 - Numerical Analysis 



Unclassified 

- Unlimited 


19. Security Classif. (of this report) 

Unclassified 

20. Security Classif. (of this page) 

Unclassified 

21. No. of pages 
22 

22. Price 

A0 3 


NASA FORM 1626 OCT 86 


NASA-Langley, 1991 




